! Huistaak 3 TWS
!
! Roel Matthysen
! s0202264
! 1e Master Wiskundige Ingenieurstechnieken

module mod_calc_roots

	implicit none
	public
	
	! Parameters voor enkele en dubbele precisie
	integer,parameter :: sp = kind(0e0)
	integer,parameter :: dp = kind(0d0)
	
	contains

	subroutine calc_roots_single(b,c)
		real(sp),intent(in) :: b,c
		real(sp) :: x1(2),x2(2),D
		
		D = sqrt((b/2)**2-c)
		x1(1) = -b/2+D
		x1(2) = sign(abs(b/2)+D,-b)
		x2(1) = -b/2-D
		x2(2) = c/x1(2)

		write(*,*) '-------------'
		write(*,*) 'Single Precision'
		write(*,*) '-------------'
		write(*,'("D = ",es13.6)'),D
		write(*,*) 'Algoritme 1'		
		write(*,'(4x,"x1 = ",es13.6)'),x1(1)
		write(*,'(4x,"x2 = ",es13.6)'),x2(1)
		write(*,*) 'Algoritme 2'		
		write(*,'(4x,"x1 = ",es13.6)'),x1(2)
		write(*,'(4x,"x2 = ",es13.6)'),x2(2)

		
		
	end subroutine
	
	subroutine calc_roots_double(b,c)
		real(dp),intent(in) :: b,c
		real(dp) :: x1(2),x2(2),D
		
		D = sqrt((b/2)**2-c)
		x1(1) = -b/2+D
		x1(2) = sign(abs(b/2)+D,-b)
		x2(1) = -b/2-D
		x2(2) = c/x1(2)
		write(*,*) '-------------'
		write(*,*) 'Double Precision'
		write(*,*) '-------------'
		write(*,'("D = ",es22.15)'),D
		write(*,*) 'Algoritme 1'
		write(*,'(4x,"x1 = ",es22.15)'),x1(1)
		write(*,'(4x,"x2 = ",es22.15)'),x2(1)
		write(*,*) 'Algoritme 2'		
		write(*,'(4x,"x1 = ",es22.15)'),x1(2)
		write(*,'(4x,"x2 = ",es22.15)'),x2(2)
		
	end subroutine
	
end module


program quadratic_equation
	use mod_calc_roots
	
	implicit none
	! Verschilllende precisie b en c
	real(dp) :: b,c
	real(sp) :: b_sp,c_sp
	
	! De verschillende decimale precisie afdrukken
	print *,''
	write(*,'("Aantal beduidende cijfers enkele precisie =  ",i0)'),precision(1.0_sp)
	write(*,'("Aantal beduidende cijfers dubbele precisie = ",i0)'),precision(1.0_dp)
	print *,''
	
	! b en c inlezen
	print *,'Geef de waarde voor b in : '
	read *,b
	print *,'Geef de waarde voor c in : '
	read *,c
	
	! Initialisatie
	b_sp = real(b,kind=sp)
	c_sp = real(c,kind=sp)
	
	! Wortels berekenen in enkele precisie
	call calc_roots_single(b_sp,c_sp)
	! Wortels berekenen in dubbele precisie
	call calc_roots_double(b,c)
	! -------------	
	! BESLUIT
	! -------------	

	! Bij het uittesten van b=-577 en c=13 werd de volgende uitvoer verkregen:

	! Aantal beduidende cijfers enkele precisie =  6
        ! Aantal beduidende cijfers dubbele precisie = 15
	! 
 	! Geef de waarde voor b in : 
	! -577
 	! Geef de waarde voor c in : 
	! 	13
 	! -------------
	! Single Precision
 	! -------------
	! D =      2.88477E+02
	! Algoritme 1
	!    x1 =     5.769775E+02
	!    x2 =     2.252197E-02
	! Algoritme 2
	!    x1 =     5.769775E+02
	!    x2 =     2.253121E-02
	! -------------
	! Double Precision
	! -------------
	! D = 2.884774687908919E+02
	! Algoritme 1
	!    x1 = 5.769774687908919E+02
	!    x2 = 2.253120910813777E-02
	! Algoritme 2
	!    x1 = 5.769774687908919E+02
	!    x2 = 2.253120910811763E-02

	! Afgaande op deze resultaten is het duidelijk dat de precisie op x1 in enkele precisie voor beide algoritmes gelijk is.
	! Voor Algoritme 2 is x2 gelijk aan de referentiewaarde berekend met dubbele precisie. Voor algoritme 1 wijkt deze 
	! waarde af. In deze situatie is algoritme 2 dus numeriek het meest stabiel.
	! Een uitzondering is bijvoorbeeld wanneer b = c = 0. Dan zal algoritme 2 falen omdat er een deling door 0 wordt 
	! uitgevoerd.

end program
